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Abstract 

We present a stable and convergent method for studying a system of gas and dust, 
coupled through viscous drag in both non-stiff and stiff regimes. To account for 
the effects of dust drag in the update of the fluid quantities, we employ a fluid 
description of the dust component and study the modified gas-dust hyperbolic sys- 
tem following the approach in Miniati & Colella (2007). In addition to two entropy 
waves for the gas and dust components, respectively, the extended system includes 
three waves driven partially by gas pressure and partially by dust drift, which, in 
the limit of vanishing coupling, tend to the two original acoustic waves and the 
unhindered dust streaming. Based on this analysis we formulate a predictor step 
providing first order accurate reconstruction of the time-averaged state variables at 
cell interfaces, whence a second order accurate estimate of the conservative fluxes 
can be obtained through a suitable linearized Riemann solver. The final source term 
update is carried out using a one-step, second order accurate, L-stable, predictor 
corrector asymptotic method (the a-QSS method suggested by Mott et. al. 2000). 
This procedure completely defines a two-fluid method for gas-dust system. Using the 
updated fluid solution allows us to then advance the individual particle solutions, 
including self-consistently the time evolution of the gas velocity in the estimate of 
the drag force. This is done with a suitable particle scheme also based on the a-QSS 
method. A set of benchmark problems shows that our method is stable and conver- 
gent. When dust is modeled as a fluid (two-fluid) second order accuracy is achieved 
in both stiff and non-stiff regimes, whereas when dust is modeled with particles 
(hybrid) second order is achieved in the non-stiff regime and first order otherwise. 

Key words: Godunov methods, Particle-In-Cell methods. Stiff equations 
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1 Introduction 



We wish to solve the system of partial differential equations describing two 
systems coupled through the exchange of momentum through a viscous term 
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proportional to their relative velocity. This situation characterizes a variety of 
problems, among others gas-dust coupling in protoplanetary disks, the motion 
of polymer molecules in biological fluids [26] , drift of different ion species in 
the planetary plasma [5]. We are interested in addressing the case in which 
the viscous coupling becomes stiff, such that the relaxation time character- 
izing it is significantly shorter than the smallest timescale characterizing the 
fluid system, typically defined as soimd crossing time of a resolution element. 
Although the results presented in this paper can possibly be extended to other 
physical systems as those mentioned above, in the following we specialized our 
analysis to the case of a gaseous and a dust component, coupled through a 
drag term as well as gravity, typical of a protoplanetary disk. Without loss of 
generality we start considering the problem in one dimension. The gas com- 
ponent is described by the equations of hydrodynamics with a suitable source 
term, namely 
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The fluid quantities (with subscript g for gas) have their usual meaning, 
whereas fd and — V0 describe drag and gravitational acceleration, respectively, 
which will be specified below. The dust particles move along the following tra- 
jectories in phase-space 



^=Vd, (3) 

dvd , , „ , , . , 

— = -Kd [Vd - Ug) - V0, (4) 

where, Xd and Vd are the dust particles position and velocity, respectively. We 
consider particles sizes that are small compared to the gas particle mean free 
path (Epstein's regime) so that the drag coefficient is 

i^d = nopgC, Kq = ^. (5) 

PdS 

Here, pd and s indicate the dust grain's mass density and size, respectively, 
and c is the gas speed of sound. As shown below (Sec. 3.1), the back-reaction 
exerted by the dust particles on the gas in Eq. (2), takes the form 

fd = -K.g{Ug - Ud), Kg = HoPdC, (6) 
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where Ud is the average dust velocity in the neighborhood of the considered 
fluid element, also to be defined below. Eq. (5) and (6) indicate that the 
dust-gas coupling coefficients become very large in the limit of either/both 
small grain sizes or/and large dust density. In this case the relaxation between 
gas and dust is fast compared to the sound crossing time, which leads to 
numerically stiff conditions. In addition, in the case of large dust densities, 
the dust back-reaction on the gas dynamics is considerable and needs to be 
accounted for accurately. 

In addition to the stiff conditions mentioned above, the problem is compli- 
cated by the fact that on the one hand the dust particles are collisionless 
and therefore best described with a particle method in phase-space. However, 
in the stiff regime, the particles effectively interact on a short timescale not 
only with the surrounding gas but also with the surrounding particles. Taking 
into account such interactions using directly a particle description complicates 
considerably the problem, to a level that in fact is not tractable. 

In this paper, we formulate an algorithm to address such a problem. In partic- 
ular, for the purpose of modeling efficiently the effects of the drag force by the 
dust particles on the fluid, we flrst obtain fluid description of the dust particle 
component, using a Particle-Mesh method. We then define an extended con- 
servative system that includes both the gas and dust variables. This systems 
is advanced one time step following a method which is an extension of the one 
described in Miniati and Colella [19]. More specifically, we use the fiuid de- 
scription of the dust to account for the modifications of the drag terms on the 
hyperbolic structure of the gas equations. This allows us to formulate a pre- 
dictor step that gives first order accurate reconstruction of the time-averaged 
state variables at cell interfaces, whence a second order accurate estimates of 
the conservative fiuxes can be obtained. With the dust component still de- 
scribed as a fluid, and a second order estimate of the fluxes for our extended 
conservative system, the source update is flnally carried out using a one-step, 
second order accurate, predictor corrector asymptotic method as suggested by 
Mott, Oren and van Leer [21]. At the end of this procedure we have obtained a 
fiuid solution that includes self-consistently the collective effect of the particles 
drag on the fluid. 

We can then update the individual particle solutions using a suitable parti- 
cle scheme that follows the particles along their characteristic trajectories in 
phase space, taking into account the effect of drag from the fluid. The above 
knowledge of the fluid solution at the current and next time-step is crucial, 
however, to include self-consistently the time evolution of the gas velocity in 
the estimate of the drag force. 

Finally, dust and gas are also coupled through gravity. However, this can be 
achieved in the standard way (e.g. [19]), namely by applying to both dust and 
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gas components the gradient of their cumulative potential, 0, defined by the 
following Poisson's equation 



Several codes have already been published in the literature in order to study 
the coupled dust-gas-dynamics in protoplanetary disks. Among others, nu- 
merical codes based on the Smoothed Particle Hydrodynamics method (SPH) 
approach have been developed for multifluid system [20,3,16], in addition to 
the simpler case in which dust is treated as a test particle and backreac- 
tion is neglected [24]. Grid based two fluid methods have been around for a 
long time [8]. More recently they have been extended to the MHD case [11], 
and to one-fluid models in which dust and gas are perfectly coupled [2]. Grid 
based hybrid (fluid+particles) methods have also been developed, including 
higher oder (sixth) spectral methods [14] and higher order (second) Godunov's 
method [1] . The novelty of our method consists in its ability to handle a vari- 
ety of numerically stiff conditions both in the two-fluid and hybrid approaches. 
This is relevant because in a realistic setting stiff conditions arise in limited 
regions of a system and at a certain point in time during its evolution. 

This paper is organized as follows. In Sec. 2. the particle integration scheme 
is presented. In Sec. 3 we describe the two-fluid approach that allows us to 
update the fluid solution one time step. This section, therefore includes a 
description of the methods for the fluid treatment of the dust component, 
the details of the semi-implicit predictor-corrector used for the source update. 
In Sec. 4 we present the Godunov predictor step, derive the characteristic 
analysis and provide a linearized Riemann solver, for the extended gas-dust 
fluid system. The stability of the method is addressed in Sec. 5. Accuracy and 
convergence tests are presented in Sec. 6 and a summary with discussion in 
Sec. 7 concludes the paper. 



2 Pctrticle Integration Scheme 

The particle positions and velocities is updated with a predictor corrector 
method based on a variant of the Quasi-Steady-State method proposed by 
Mott, Oran and van Leer (the so called a— QSS, [21]). Because we use this 
scheme extensively in this paper, we first describe it briefiy in the following. 

Consider the first order ordinary differential equation (ODE) 



(7) 



dy 
dt 



■p{t,y)y + q{t,y), 



(8) 



y{to) = 2/0, 



(9) 
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with 1/ e R, p, g : M X M ^ M. If p and q are constants, the above system has 
the following exact solution given by Duhamel's formula, 

?/(i) = yoe-^* + ^(l-e-^'*). (10) 

QSS methods are based on the asymptotic behavior of the above solution. 
When p and q depend on (y,t), a first-order algorithm can be obtained by 
setting p = po = p{to), q = qo = q{to). This approach corresponds to the 
simplest QSS method and can be used to develop higher order methods by 
incorporating to some degree the time dependence oip,q into the solution. All 
QSS methods, however, reproduce the exact solution as p, q become constant. 
There are several QSS derived methods in the literature (cf. [23] for a review). 
Here we shall adopt the a— QSS method of Mott et al. [21]. This is a single 
step, second order accurate, A-stable predictor-corrector method that can be 
summarized with the following procedure: 

.fa + Ai)^.o + At ^/;-^*^^ . (12) 

q* = qo + a{pAt) [q{y, to + At) - go], (13) 

Eq. (11) and (12) correspond to the predictor and corrector step respectively. 
In addition to the case where p, q are constants, the above algorithm also 
returns the exact solution when p is constant and q is linear in time or p is 
linear in time and g = [21]. 

In the following we derive a predictor-corrector algorithm based on Eq. (11)- 
(12) to integrate the equation of motion of the dust particles given by (3) and 
(4). In the absence of drag it effectively reduces to a kick-drift-kick variant 
of the leapfrog scheme. For its derivation we assume that the gas velocity 
solution, Vg, has already been computed by the fluid scheme described below, 
at both and t""*"^ = t" + At. The drag coefficient in (4) is a function of the 
gas density and sound speed and is interpolated at the particle position with 
a Particle-Mesh method, i.e. 

i^di^p) = ^w[xp - x{i)]Kd{i), (15) 
p 

where w{s) is a weight function, x{i), is the i-th cell's center, and the summa- 
tion is carried out over all cells. The velocity of the gas entering Eq. (4) must 
also be obtained by way of interpolation. By using the following expression 

'"^^^^^ ^ ^^^50 ^^^^^ ~ ^(i)]'«<i(i)«fl(i)> (16) 
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total dust+gas momentum is conserved by construction in the non-stiff limit, 
and simple tests show that it is also conserved with high accuracy in the stiff 
regime. 

We now present the algorithm, but a more detailed description of its derivation 
is provided in Appendix A. We first predict the particle position at t"'~^^ as 

Xd = x'^ + (3{K'^At)v2At + [1 - (3{K-,At)]vl^At + ^V^ At^, (17) 

where = ^—j—- As the particle travels from x^ to Xd the gas velocity 
changes from Vg — Vg{xd) to Vg'^^ — Vg{xd)- Since the above change is partially 
due to drag relaxation and partially due to the motion of the particles across 
a velocity gradient, we assume that the velocity experienced by the particle 
evolves as follows: 

vg{t) = < + A^, (e-'^^^*^ + 1 - e-'^^*) , (18) 

where, Av„ = f "^^ — f The above expression allows for a relaxation of the 
gas velocity towards the value v^^-^ much faster then linear, which is important 
in the stiff case {KgAt ^ 1). For the stiff case, this introduces some differences 
in our integration scheme with respect to the pure a-QSS method. 

With the assumption of Eq. (18) we can carry out the time integration of 
the equation of motions for the particle velocity and position. In doing so, as 
in the a-QSS approach, we replace the time dependent drag coefficients with 
their time averages, that is the average between the values at position Xd and 
Xd, Ks ^ Kg — [K's{xd) + K'si^d)]/'^, s — d,g. We thus obtain 

^T' = + v2^t + [1 - /3 (KdAt)] (< - v2) At + 



— KgAt 



+Avg { e-'^o 



1 1-p (Kd^t) 



2 

^ ^ Kdf^jKg^t) - KgPjKd^t) \ ^^g^ 



2 KdAt 

v2^' = v2 + P{>id^t)Kd{v^ - v2)At 
+Avg |e-'^"«^*[l - (3{Kd)] + -^^[PiKgAt) - P{KdAt)]At] . (20) 

I Kd-Kg J 



Remarkably the above particle method contains no explicit term arising from 
the stiff coupling of the particle component. Of course, the /3 terms take into 
account that the gas-dust coupling is fast compared to the timestep At. How- 
ever, each particle motion is integrated individually though it effectively de- 
pends on the other particles solutions, and the scheme is essentially explicit 
in time, as it only involves the particle solution at time t = nAt. Still the 
scheme is stable and convergent, and this is due to the fact that the gas veloc- 
ity solution at times t and t + At, with which the particles interact, already 
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contains the effect of the dust component to second order accuracy, even in 
the stiff regime. The gas velocity solution is in fact the only quantity entering 
Eq. (19)-(20). We have actually tried to employ more sophisticated approaches 
than Eq. (18), which would also involve the dust fluid velocity, but despite the 
higher degree of complexity and computational cost, they did not improve on 
the algorithm accuracy. 



3 Two-Fluid Semi-Implicit Predictor-Corrector 



3.1 Two-Fluid Description 



In order to efficiently model the collective effect of the drag force exerted 
by the dust particles on the fluid we use a fluid like description of the dust 
component. Thus, using a Particle-Mesh method we define the dust density 
and velocity field on a Cartesian grid as 



Ud{i) = ^w[xp - x{\)]vp, (22) 
p 

where as in the previous section w{s) is a weight function, a;(i), is the i-th 
cell's center, but now the summation is carried out over all the dust parti- 
cles. Note that the position and velocity of the dust particles are updated 
independently at each time-step using a particle-method described in Sec. 2. 
The time evolution of and PdUd is then obtained from the zero-th and first 
velocity moments, respectively, of Boltzmann's equation for the dust distri- 
bution function, gd{x,v,t). Using the particle equation of motion (3)- (4), and 
the appropriate collision term we obtain 



|^ + |te»J = 0, (23) 
d d 

^ (PdUd) + (PdUdUd) = -pdK'diud - Ug) - pdV(t>, (24) 

where closure is granted by neglecting higher order velocity terms as suggested 
by [13]. However, different closures can in principle be employed starting from 
the particle description as need be. Finally, by Newton's third Law, Eq. (24) 
indicates that the collective drag force exerted on the gas by the dust particles 
ought to be 

fd^ -Hg{Ug -Ud), Kg^Kd — , (25) 

P9 
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ensuring momentum conservation. 

We can now extend the set of conservative variables, U, to include the density 
and momentum of the dust component, 



U = (P, pUg, PE)'^ {Pg, PgUg, pE, ^rf , Prft^d)^- 



(26) 



The extended set of equations for U is obtained by combining (l)-(2) and 
(23)-(24). After rearranging the source terms, it reads 



(27) 



where 



F([/) = 



pUg 

pul + P 

Ug [PE + P] 

pUd 





UgPgW(f) 




^0 ^ 

—Kg Kd 







-Kd J 



(28) 



3.2 Semi-Implicit Method 



Given our extended system of equations (27) we aim for a scheme in which 
an explicit approach is retained for the non-stiff conservative hydrodynamic 
term, V ■ F, and a semi-implicit method is employed for the stiff part of the 
source terms. As in Miniati & Colella [19], our time discretization for the 
source terms is a single-step, second-order accurate scheme. However, instead 
of a method based on the deferred correction approach [9], here we derive our 
predictor corrector using again the «— QSS [21]. The reason is that although 
stable and convergent, unlike the «— QSS method, the deferred correction 
method is not L-stable (although a combination of two such methods can 
lead to L-stability) . This can lead to large errors, particularly with regard to 
momentum conservation, in a hybrid method in which fluid and particles are 
stiffly coupled. 

Our semi-implicit method consists in solving the following collection of ODEs, 
one at each grid point, 

^ = KuU + S{U)-{V-Fr+'^, (29) 
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where we view the time-centered flux divergence as a constant source, whose 
computation using a modified Godunov method is described below. Our predictor- 
corrector step then reads: 



U{to + At) = e^»^*f/o + X^„(At) [SiUo) - (V ■ F)"+^ 
U{to + At) = e^^'Uo + I'kiAt) [S{Uo) - (V ■ F)"+^ 

S{U) - S{Uo) At^ 



At, 
At 



+Xj,(At)- 



At 



(30) 



(31) 



where, K = [K{Uo) + K{U)]/2, is used in Eq. (31), and we have defined the 
set of operators 



(32) 



Note that gravity is actually unaffected by the operators X3(t); however the 
form of Eq. (30)-(31) is suitable for a more general form of S{U), which may 
even include stiff terms associated, for example, with an endothermic source] ^ I 
Although the form of predictor corrector in the above equations appears differ- 
ent from the original system (11)-(12), their equivalence can be easily verified 
by replacing the operators K, Q with the coefficients p, q. It is also easy to see 
that in the non-stiff limit, the above scheme reduces to the usual second order 
accurate explicit formulation 



U{to + At) = (1 + KAt)Uo - At (V ■ F)"+^ + ^ \s{U) + 5(f/, 



(35) 



4 Two-Fluid Predictor Step 



In order to compute our predictor step, we cast the extended two-fluid gas-dust 
system in primitive form as follows: 

dW dW 

-^ + AiW)^ = KW + G, (36) 



^ In this case, if S indicates the endothermic source, Eq. (29), as well as Eq. (30)- 
(31), should be modified by replacing 

S{U)^S{U) + J:{Uo), (33) 
K{U) ^ K{U) + Vu^iU). (34) 
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where the primitive variables are 



^ = {Pg,Ug,P,Pd,Ud)'^ , 



(37) 



and the quasi-hncar operator, ^4(1^) = V[/VF • V[/F • VwU, which includes 
an advection and a Lagrangian component, is 



A{W) = Ugl + AL, Al^ 




p-' 

Q Pg(? 

{ud-Ug) 






Pd 



\ 



1^0 {ua-Ug)] 



(38) 



Finally 



K 





—Kg Kg 







G 



( ^ 






(39) 



yo -Ka j 

represent the drag operator and the residual source term, including gravity, in 
primitive form. Following the method in Miniati & Colella [19] we can follow 
the dynamics of the system along Lagrangian trajectories 



DW dW 

-Al{W)^ + KW + G, 



Dt ' dx 

with solution given by Duhamel's formula: 



(40) 



W{t) 



;^*iy(0) - /* 



.K{t-T) 



OX 



dr. 



(41) 



This allows us to define a modified dynamics reading 

DW dW 

+ Ij,{t)AL^ = lK{t) [KW{Q) + G] + 0{t). 



Dt 



dx 



(42) 



4-1 Characteristic Analysis 



We use the quasi- linear system (42) with t = At/2 to compute the Godunov 
predictor step. Thus we analyze the modified hyperbolic structure of that 
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system of equations. With the choice of K in (39), from Eq. (32) we obtain 

/i n n n n \ 



1 

K K 

10 

1 



V' 

where k — Kg -\- K^is the total relaxation rate, and 



(43) 



where 



/l 
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Kg(l-/3) 
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K 
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1 




o<p< 




i/tAt 





(44) 



(45) 



We can then define the modified linear operator 



lK{At/2)AL = 



eff 



P, 







KPg 

>2 






















{Ud - Ug) Pd 

^-^{Ud-Ug)j 



(46) 



with associated characteristic equation 

.^Kg^ ^Kd 



A [A - {Ud - Ug)] 



K 



{Ud -Ug) - X 



l^d + ^l^g 2 , n 2/ \ 
''C + /3c [Ud-Ug) 



K 



0. 



7) 

The system eigenvalues are given by the roots of the above equation. The 
first obvious solution, A° = 0, correspond to an entropy wave for the gas 
component. Since an entropy wave does not involve velocity perturbations, 
we expect this type of wave to be unaffected by the presence of a drag term. 
The other obvious eigenvalue, A'^ = [ud — Ug), is the equivalent of an entropy 
wave, but for the dust component. Perturbations in propagate only along 
the characteristic curve associated to this eigenvalue and, therefore, do not 



11 



affect otlicr primitive variables. For tfiis reason, in tfie following we simplify 
the characteristic analysis by dropping out the components and assume 
that: pd is transported as a passive scalar with speed A'^ and its intermediate 
state Riemann solution is simply given by the average of left and right state 
values. 

The remaining eigenvalues of the system are given by the roots of the cubic 
polynomial appearing in Eq. (47), which can be shown to be always distinct 
and real, with a few exceptions discussed below. The set of relevant eigenvalues 
is then 



A° = 0, (48) 

(49) 













73 


71 


( Sj^ 



(f 2 
3 + 3' 



(p 4 
1+3^ 



(50) 
(51) 



where 



C, OU = — dU, OU — Ud — Ug, 



K 



K 



if = COS" 



5u 



25u + 9c2 (l - 



(52) 
(53) 



and cos^^ conventionally indicates the principal value of the multivalued in- 
verse cosine function. With this convention we always have A+ > (A^,A°) > 
A~, whereas all three cases A'^ ^ A^ are admitted. An extra eigen-mode 
driven by the dust drift has appeared. Under stiff dust-gas coupling condi- 
tions (/3 < 1), this mode mixes with the pressure driven modes thus acquiring 
an acoustic-like character. Likewise, the properties of the original sound waves 
are contaminated by the dust component. The degree of mutual contamina- 
tion is regulated by the 'mixing angle' (p. In general sound speed is affected 
by both the additional inertia and by the motion of the dust particles. To see 
this in more details in the following we explore the asymptotic behavior of the 
waves in a few cases of interest. The general results are summarized in Fig. 1. 
In the non stiff limit, we have 



lim A"^ = c. 



lim A = 



lim A^ = 5u, 



(54) 
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Fig. 1. Solid lines: eigenvalues versus {ug — Ud), in sound speed units, for /3 = 0.25 
and Kg/ K = 0.65. The dash lines indicate the various asymptotic branches corre- 
sponding to A = c (high, horizontal) A = — c (low horizontal) and \ = 6u (diagonal). 



thus recovering the expected values. However, as the dust particle drift van- 
ishes we have 

(3k, 

lim A"*" = c, lim A~ = — c, lim = — 5u. (55) 

Su^O Su^O Su^O Kfi + pKg 



Here c is given in Eq.(52) and correspond to a sound speed in which the 
gas density is replaced by an effective average between the dust and the gas 
densities. So, in the stiff limit (/3 — )• 0), it is the total density that enters the 
definition of the sound speed (which can be dominated by the dust density). 
Similarly, in the limit of large values of Su we obtain (see dashed lines in Fig. 1) 



lim A^ 

Su—^oo 

lim A^ 

5u— >— oo 



lim A 

5u—^oo 

lim A~ 

5u—^ — OD 



6u, 



lim A^ 

5u—^oo 

lim A^ 



-c. 



(56) 
(57) 



where we have used 



\ Kg + f3Kd 



(5J 



which represent an effective sound speed when the perturbation propagates 
opposite to the direction of the dust particles velocity. The said propagation 
speed vanishes in the stiff limit, but it always coincide with the sound speed 
(c) in the limit of negligible dust density (i.e. Kg — )■ 0, Kd ^ k). 



To conclude the analysis we introduce the array of right eigenvectors (by 
column) 
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R = 



^ ^ 

Pg Pg 





pg 
„2 



{c~\-){c+\-) (£-AX)(c+A>') (£-A+)(c+A+) 



(59) 



where, to simplify the notation, we have introduced the reduced density, /j, = 
PgPd/{Pg + Pd)- Similarly, the array of left eigenvectors (by row) is 



(A++Ax)p, 


c2+A+AX 


H5u{l3-1) 


(A+-A-)(AX-A-) 


c2(A+-A-)(AX-A-) 


(A+-A-)(AX-A-) 





1 





(A++A")p3 


c2+A+A- 




(A+-A=<)(AX-A-) 

(a>^+a-)p3 


c2(A+-Ax)(Ax-A-) 
c^+A^A- 


{\+-\x)(XX-\-) 
^i{ua-Ug){l3-l) 



(60) 



(A+-A-)(A+-Ax) c2(A+-A-)(A+-A>' 



Just like for the eigenvalues, it can be easily verified that the definitions (59) 
and (60) tend to the usual expression for the left and right eigenvectors in the 
non stiff limit. 



4.1.1 Loss of Strict Hyperbolicity 

When either one of the following three cases occurs, /3 — )■ 1, 5m — )■ 0, — )■ 0, 
(implying either pg = or pd = 0) the first, third and fourth right eigenvectors, 
appear to become singular. However, the asymptotic analysis shows that in 
these limits the expressions, =F c, approach zero quadratically in (/3 — 1) 
and linearly in both pa and 5u, so that the first and fourth right eigenvectors 
are always well defined. 

While the third right eigenvector still diverges, the corresponding left eigenvec- 
tor tends to zero at the same rate (because A+-|-A~ is quadratic in (/3— 1) and 
linear in both pa and Su, whereas + A'^A", is quadratic in both (^ — 1) and 
6u, and linear in p^), so that the characteristic decomposition is non-singular 
in the above limits. Note that this "illness" of the third right eigenvector could 
be formally cured by renormalizing the third left and right eigenvectors by the 
factor fs and f^^, respectively, where 

- /.5m(/3 - 1) ■ ^^^^ 



The eigensystem can still become singular when two eigenvalues become iden- 
tical, i.e. strict hyperbolicity is lost. Inspection of Eq. (49)-(51) and (53) in- 
dicates that it is possible to have A~ = A^ or A+ = A^ when (f) — or (f) — 
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respectively. These relations are satisfied when 5u = ^c, respectively, and si- 
multaneously either /3 = 1, or = is verified. As it appears from (59) and 
(60), in this case the third (x) eigenvectors (left and right) become parallel 
to either the first (— ) or fourth (+) eigenvectors, respectively. When this is 
the case we modify our characteristic synthesis as suggested by Bell et al. [4] . 
For example, if ak = J2k h ■ SW are the coefficients of the characteristic de- 
composition of a perturbation 6W with respect to the set of left eigenvalues 
Ik, when — < e, with typically, e ~ 10~^, the above perturbation would 
be reconstructed as 

SW akVk ay,r± + ^ akVk, (62) 

k k^x 

where rfc is the set of right eigenvectors. 



4-1.2 Addition of Endothermic Processes 

We can easily generalize the above results to include stiff endothermic pro- 
cesses. Suppose the rate of change of the gas internal energy, e = P/ Pg{'j — 1), 
is 

de 

^ = A(e,p,). (63) 
In this case the modified linear operator defined in Eq. (46) becomes 






apgC^ 











KPg 

(l-a)ApPg Q 



Ae 








^ 

^-^{ua-ug) 


{Ud - Ug) Pd 

^-^{Ud-Ug)l 



Kd(l-/9) 

KPg 



(64) 



where 



dk , OA 

Ae = ^, = — 

oe opg 



- 1 



A^At 



-, < a < 1. 



(65) 



Then it can be shown that the above characteristic analysis remains valid 
provided the gas sound speed is replaced with the following effective value [19] : 



Ceff 



l + a(7-l)-(l-a) 



KPg 

AeS 



P9. 



(66) 
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Fig. 2. Four-waves structure of the Riemann problem. 

As discussed at length in [19] the above scheme is apphcable as long as the 
fluid is thermally stable [10], namely: 



> 



a 



ApPg a(7-l) + l' 



(67) 



4-2 Linearized Riemann Solver 



The characteristic analysis derived in 4.1 defines the structure of the solu- 
tion to the Riemann problem for the gas-dust system. Unlike ordinary hy- 
drodynamics, this is now characterized by three acoustic wave families and 
one entropy wave propagating at the fiuid velocitjIZl- The corresponding four 
characteristic curves determine in general five distinct regions, corresponding 
to the initial left (L) and right (R) states and three intermediate states of 
the gas. Fig. 2 illustrates this, for the case in which < A*^. In the x — t 
plane, the solution cone is separated by the unperturbed R and L states by 
the fastest and slowest eigenvalues, respectively. The R-state is connected to a 
*R-state by a rarefaction fan or a shock depending on whether Aj^ is faster or 
slower than A^^. The same applies to the connection between the L-state and 
*L-state. Then if A^ < A" as assumed in Fig. 2, the *R-state is connected to 
a #-state by an ordinary contact discontinuity, where the density jumps but 
pressure and velocity remain constant. Finally, another acoustic wave sepa- 
rates the #-state and *L-state, which again either takes the form of a shock 
or a rarefaction wave depending on whether the characteristics with speed A^^ 
and A^ converge or diverge. 

In order to calculate the intermediate states described above, we use a lin- 

^ For the sake of clarity, in the following we continue to ignore the contact discon- 
tinuity in the dust component, which travels at speed A'^ without interacting with 
the rest of the primitive variables. 
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earized Riemann solver [27], which is discussed in detail below. Then, in gen- 
eral, given the input left and right states Wl,Wr, we can use the left and 
right eigenvectors, {lk,'>^k}k=i...n and the corresponding eigenvalues A''', com- 
puted at some intermediate state, to decompose the perturbation along the 
characteristics as follows: 



k 

E 

k'=l 



Then the solution to the Riemann problem along the ray x/t — 0, Wrp, is 
given by 

\k 

Wrp ^Wl+ Y1 ^krk + ^kTk + TT^Tk^krk, 

max(A|,A|)<0 A* >0>A^,A|+A^<0 A*<0<A^ ^ ^ 

The first sum is the correction to the state for the waves that are unam- 
biguously have a negative speed; the second is a correction to the state for 
transonic shocks; and the third a correction for transonic rarefactions, that 
prevents the formation of entropy-violating shocks. In the following we use 
a slight variation of this approach that is better able to preserve reflection 
symmetries. 

Recall that we have three acoustic-like wave families A; = +, — , x with speeds 
A*^ given in (49)- (51) and the remaining wave propagates at the fluid speed, 
A°. Assume for the moment that we can compute the eigenvalues A'^ at the 
intermediate states, *L, *R, following a procedure described in detail be- 
low. Then the approximate Riemann solver is given as follows. 

If A° > 0, then 

if min(A^^, A^) > or A^^ > > A^ and A^^ + A^ > then 



Wrp^ 



'W^L if max(Ai, A,i) < or A^ > > A,^, A^ + A,^ < 
Wl + T^{W,L - Wl) if AZ < < a;^ 
Wl otherwise 



else if A^^ < < A| then 

Wrp = W.L + Tf^{W# - W.l) 
else 



Wrp = W^ 
endif 

else if A° < 0, then 
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if max(A|, A^^) < or > > A^^ and A^ + A^^ < then 

r W,R if max(A+j, A+) > or A+, > > A+ A+, + A+ > 
Wrp ^\wr + -j^iW^R - Wr) if A+^ < < A+ 
\Wr otherwise 

else if < < A^^ then 

Wrp = W^R + T3%r(W^# - W,r) 

else 

Wrp = W# 

endif 

endif 

In order to complete our description of the approximate Riemann solver, we 
now outline the method for computing the intermediate states. Basically we in- 
tegrate the characteristic equations along the characteristic directions to relate 
the jumps experienced by the primitive variables across adjacent states [12]. 
The characteristic directions are assumed to be constant and are computed at 
the foot of the characteristics. The characteristic equations are obtained by 
setting ■ dW — for each eigenstate and they read 



dp - (A^ A+)rfM, - X, (A^ A+)(iMd = 0, (68) 

dp — c^dpg = 0, (69) 

dp - (A-, X+)dug - Xd (A", A+)rfn, = 0, (70) 

dp - x+(A^ A-)rfM, - xl(A^ A-)rfMd = 0, (71) 

where the characteristic slopes in phase-space are given by 



x.^(A^A-)=^^^^p.c^ (72) 

Xd(A ,A )-^^^^5-^c, (73) 

and the analogous expressions for x^di Xgd obtained by permutation of the 
symbols +, — , x in the above equations. The above procedure of integrating 
the characteristic equations yields 12 equations in 12 variables. However, nei- 
ther the gas or dust fluid velocities nor the gas pressure change across the con- 
tact discontinuity, reducing the system to 9 equations in the following 9 vari- 
ables: Pg,*Lj '"g,*Lj Pg,*L, '"(i,*Lj Pg,*R, Ug^^R, Pg,*R, U^^^r, Pg^#, with Pg,#, 

coinciding with the corresponding variable of either the *L or *R state. In order 
to keep the system linear we make two further approximations when comput- 
ing the characteristic directions (72)-(73) in the intermediate states: (1) we 
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assume that the gas-dust velocity difference, Ug—Ud, in the intermediate *L,*R- 
states, can be approximated with the values corresponding to the L.R-states, 
respectively. And (2), that, when connecting the #-state with the *L-state 
(*R-state) via the — (+) characteristic, i.e. when < A° (A^ > A°), we can 
use the density of the *L-state (*R-state) instead of the value corresponding 
to the #-state. With this latter approximation the relations between velocity 
jumps across the characteristic line x, remain unchanged whether A^ < A° or 
A^ > A°. The solution then reads 



Ug,*R = Ug^R + 



P9,#- 



{Pl - Pr){1 - ^) + Kl - Ug,n)ieL + ^Vl) + {ug^L - ^d,RY-^ 



P*R—Pr + dR{Ug,R — Ug,*R), 

Ud,*R^Ud,R + r)R{'^g,R - 

_ , P*R — Pr 
Pg,*R - Pg,R H :o , 



(75) 
(76) 

(77) 



Ug,*L = Ug,L + 



{Pl - Pr){1 + {ug,L - Ug,R){0R + '^m) + {ud,L - ^'i,RY-^ 

{eL-9R)-m'-i!^ + rj,'-^ 

(79) 
(80) 



P*L ^Pl + dL{Ug,L - Ug,*L), 
Ug,*L = Ug^L + r]L{Ug,L - Ug,^L): 

, P*L — Pl 

Pg,*L - Pg,L H "2 , 



^Pg,*R+^^^^ ifA,\> A^ 



P*R—P*L 



_ _ _ Pd,L + Pd,R 
Pd,*L — Pd,*R — Pd,# — 7^ > 



(81) 
(82) 
(83) 



where we have defined the following coefficients. 



^R 



Xg,RXd,R Xg,RXd,R 



Xd,R Xd,R 



A+ 



Xq,LXd,L Xq,LXd,L 



Xg,R Xg,R 
Xd,R ~ Xd,R 



p^jA+f-c^ 
A+ pSu{l3 - 1) 



R Xd,L Xd,L 

Xg,L ~ Xg,L _ 



r' Xd,L-Xd,L [A pSu{P-l)\^' 



u- 

(84) 

Pg {x-y-c^- 



Xg.*RXd,*L Xg,*LXd,*R 
Xd,*L Xd,*R 



Xg,*L Xy.*R 

Xd,*L Xd,*R 



(85) 
(86) 



The extra symbols L,*L,*R,R in the subscript of the x functions defined 
in (72)-(73) indicate the state at which the characteristic directions is de- 
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fined. Note that while the definition of the coefficients OL,Oji,r)L-,flR is based 
on the known L.R states, computing 9^, and 1]^.. involves knowledge of the in- 
termediate states, *L or *R. It was in order to avoid the nonlinearities arising 
from such dependencies that we have introduced the above approximations 
(1) and (2). As a sanity check we may note that the non-stiff limit for the 
characteristic slopes in phase space read 

limxJ = Tpc, limXd,Xp=0, limx^ = oo, (87) 
so that the coefficients in (84)- (86) tend to 

lim 6l = PlCl, fim 6r = -PrCr, hm 6^ = const., lim rjrR = 0, hm r)^ = oo, 

/3->l /3->l 13-^1 13-^1 ' 

(88) 

and the Riemann solver solutions (74)-(83) reduce to the usual expressions [12]. 



4-3 Godunov Predictor in One Dimension 



With the operator A'^^ and the sets of left and right eigenvectors that we have 
worked out in section 4.1, the Godunov predictor step is carried out as usual 
as follows. First the local slopes are defined. In particular at each point left 
and right one-sided slopes as well as cell centered slopes are evaluated and 
then a final choice on the local slope AWi is defined by using van Leer limiter. 
The upwind, time averaged left (— ) and right (-I-) states at cell interfaces due 
to fiuxes in the normal direction, q, are then reconstructed as: 



W,,^,, = ly/^ + 1 (/ - ^ ) P± ( Aiy,) , (89) 

where 

P±m= E {k-W)-rk. (90) 

±Afe>0 

The source term component is likewise accounted for as 

W,,^,, = W,,±,, + ^Xl{At/2)S. (91) 



The fiuxes at the cell faces F^^i are computed by solving the Riemann prob- 

^ n+- 

lem with left and right states given by (Wi+,Wi+i,_) to obtain W. and 

computing -Fj+i = F {w^^i'^- To modify this procedure to account for the 

effective dynamics, we use the characteristic analysis of the effective dynamics 
to perform each of the three steps. The projection operator and any limiting 
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in characteristic variables is done using the eigenvectors and eigenvalues for 
the effective dynamics derived in Sec. 4.1. Finally, the approximate Riemann 
solvers we use was defined in 4.2. 



4.4 Extension to More than One Dimension 



For directionally unsplit schemes in D dimensions an additional step is re- 
quired in order to correct the time-averaged left/right states at cell interfaces, 
Wi^±^d in Eq. (91), for the effects oi D — 1 fluxes perpendicular to the cell 
interface normal direction. Based on Eq. (42) the effect of the stiff source 
term would be accounted for by carrying out for each additional direction, g, 
a transformation 

A, ^ X°:(At/2)Ai,, + u,\ = Af, (92) 

analogous to that described in Eq. (46). In the method proposed by [7,25] the 
corrections due to transverse fluxes are computed according to a conservative 
scheme. For example in two dimension^, with q = x 

W,j,±,, = Ty,,,, - ^VuW (f^.. - Ff^ , (93) 



2Ay 



«J+2 *J 2 



where the input Wij^±^x is computed using a one- dimensional Godunov calcu- 
lation as in the previous section, the fluxes .^i. The above transfor- 
mations imply the following transverse flux corrections for the gas and dust 
velocity, 5ug, respectively: 



5ud,x 5ud,x - 



At Kg 






2Ay K 


Pg 


At Kd 


(1-/3) 




2Ay K 


Pa 



{Ud - Ug) (Uyi J 

(94) 



+1 " 



(Ud-Ug) [Uy^.J^ - Uy^,J_i) 

(95) 



where, unless explicitly indicated, all quantities are evaluated at cell center. 



^ Notation in Eq. (93) indicates that primitive variables are converted into con- 
servative variables which are then updated through conservative fluxes and then 
converted back into primitive form. 
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5 Stability Considerations 



The stability properties of the above modified Godunov's method are analo- 
gous to those discussed in [19] when considering endothermic source terms. In 
particular we note that the inspection of the characteristic analysis shows that 
the sub- characteristic condition for the characteristic speeds at equilibrium is 
always satisfied. This condition, while being necessary for the stability of our 
linearized system [27], also guarantees that the numerical solution tends to 
the solution of the equilibrium equation as the relaxation time tends to zero 



In addition, since the structure of the equations and the numerical framework, 
including the Riemann solver, remain fundamentally unaltered with respect to 
classic Godunov's schemes, we expect the usual stability conditions to apply, 
namely the familiar CFL condition on the time step 



As for the predictor corrector method described in Sec. 3.2, A-stability prop- 
erties the tt— QSS method have been discussed in detail by Mott et al. [21]. 
Although not specifically mentioned by the authors, the same stability analysis 
they present shows immediately that the a— QSS method is also L-stable. 



6 Tests 



6.1 Particle Scheme 



Before testing the convergence of the gas-dust scheme, we investigate the per- 
formance of the particle scheme presented in Sec. 2. For the purpose we con- 
sider an individual particle propagating through a medium with specified den- 
sity velocity and pressure distributions. As expected [21], our scheme repro- 
duces the exact solution when the particle moves through a uniform medium, 
and when either the coupling constant or the background fluid has linear time 
dependence. Thus in the following for our testing purposes we use a non- 
linear time dependent density, velocity and pressure field corresponding to a 
one-dimensional, right propagating sound wave solution. Note that, although 
not presented here, we obtain equivalent convergence rates as shown below 
when an external force acting on the particle is included. The wave quantities 



[6]. 





(96) 
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Table 1 

Particle Scheme Convergence Rates. 



■Insteps 


Ex 


Rx 


Ey 


Rv 


KoAt 


2 


4.6E-3 


1.4 


8.1E-4 


0.1 


5.0E1 


4 


1.7E-3 


2.2 


7.8E-4 


0.3 


2.5E1 


8 


3.8E-4 


2.2 


6.5E-4 


2.7 


1.2E1 


16 


8.4E-5 


2.0 


l.OE-4 


2.0 


6.2E0 


32 


2.1E-5 


2.0 


2.5E-5 


0.9 


3.1E0 


64 


5.2E-6 


2.0 


1.4E-5 


1.8 


1.6E0 


128 


1.3E-6 


2.0 


4.0E-6 


1.9 


7.8E-1 




% 9F,-7 


9 n 


1 nF,-fi 


9 n 


^ QP,-1 


512 


8.1E-8 


2.0 


2.6E-7 


2.0 


2.0E-1 


1024 


2.0E-8 


2.0 


6.6E-8 


2.0 


9.8E-2 


2048 


5.0E-9 


2.0 


1.7E-8 


2.0 


4.9E-2 


Table 9 

Particle 


Scheme C( )nv(n'geiic(^ 


Rates. 








■N^steps 


Ex 


Rx 


Ev 


Rv 


KoAt 


2 


4.4E-3 


1.3 


9.4E-4 


1.7 


5.0E5 


4 


1.7E-3 


2.2 


2.9E-4 


3.2 


2.5E5 


8 


3.8E-4 


2.2 


3.1E-5 


0.5 


1.2E5 


16 


8.2E-5 


2.0 


2.2E-5 


2.0 


6.2E4 


32 


2.0E-5 


2.0 


5.5E-6 


2.1 


3.1E4 


64 


5.1E-6 


2.0 


1.2E-6 


2.1 


1.6E4 


128 


1.3E-6 


2.0 


2.8E-7 


2.1 


7.8E3 


256 


3.2E-7 


2.0 


6.6E-8 


2.1 


3.9E3 


512 


8.0E-8 


2.0 


1.6E-8 


2.1 


2.0E3 


1024 


2.0E-8 


2.0 


3.6E-9 


2.3 


9.8E2 


2048 


5.0E-9 


2.0 


7.4E-10 


2.7 


4.9E2 



denoted with a 'w' subscript are given by [17] 



Uu;{x, t) = Uw,o + B sin - ^^^,0^ 
7 - lu^ 



Co + 



7 + 1 



w,0 



1 + 



2 Co 



2 

7-1 



(97) 



1 + 



7 - 
2 Co 



2t 

7-1 



(98) 
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where, cn = \ r^^^, and we use 

5 = 10-2, «^,o = 2, p^,o = l, P«,,o = l, (99) 

Using the scheme described in Sec. 2 we advance the particle position and 
velocity to a solution time tend using a progressively larger number of steps, 
A^^feps, and a correspondingly smaller time-step, At^ — tend/Nsteps- We then 
estimate the error by comparing solutions with resolution different by a factor 
2 and use Richardson extrapolation to measure the convergence rates. So, if 
s{m) is the solution obtained using m time-steps, the convergence rate is given 

by 



In 



s(m)— s(2m) 



ln(2) 

In Table 1 we report results for a choice of the coupling constant kq such 
that the propagation regime goes from mildly stiff to non-stiff. Prom left to 
right the table columns report the number of steps, the error and convergence 
rate of the particle position and velocity and, in the last column, the stiffness 
parameter, hqAI. The scheme is clearly second order accurate over the stiffness 
range reported in the Table, even though the convergence rate of the velocity 
appears to be slower for small N steps- However, this is not due to a decrease 
in convergence rate in the stiff regime. As illustrated in Table 2, the particle 
scheme is still second order accurate even for much larger values of the stiffness 
parameter. 



6.2 Convergence Rates in Smooth Flows 



In this section wc test the convergence of the method presented in this paper, 
by studying the propagation of small perturbations in a gas-dust system with a 
long and short drag relaxation timescale. We consider both a two-fluid method 
in which the dust is modeled as a fluid (and dust particles are not used), 
as well as the full fluid+particle method. In our test, the fluid component 
is initialized with uniform density, pressure and velocity values, except for 
a sinusoidal perturbation with small amplitude superposed to the x-velocity 
component. In formulae 

Ug^a; ^uo[l + Acos(27rk-r + 7r)], (101) 
P = Po = 7 = 1-4, P = Po = 0.5, Uy = Uyo = 0.7, Uq = 0.5, (102) 

where r is the position vector. Similarly the dust particles are uniformly dis- 
tributed on the grid in order to produce a uniform density, and their velocity 
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Table 3 
Run Set 



run 


D 


A 


k 


Kq 


Pd/pg 


Notet 


A 


1 


1.4 X 10-2 


(1,0) 


1 


1 


two-fluid 


B 


1 


1.4 X 10-2 


(1,0) 


10^ 


1 


two-fluid 


C 


1 


1.4 X 10-2 


(1,0) 


10^ 


10-^ 


two-fluid 


D 


1 


1.4 X 10-2 


(1,0) 


10^ 


10^ 


two-fluid 


E 


1 


1.4 X 10-2 


(1,0) 


2 X 102 


1 


two-fluid 


F 


2 


1.4 X 10-2 


(2/^/5, 1/V5) 


1 


1 


hybrid 


G 


2 


1.4 X 10-2 


(2/^/5, 1/V5) 


10^ 


1 


hybrid 


H 


2 


1.4 X 10-2 


(2/^/5, 1/V5) 


10^ 


10-3 


hybrid 


I 


2 


1.4 X 10-2 


(2/V5, 1/V5) 


10^ 


10^ 


hybrid 


L 


2 


1.4 X 10-2 


(2/V5, 1/V5) 


102 


1 


hybrid 



/i2/6nd=hydro + dust particles. 



is initialized as 

= Mo[l + ^cos(27rk-r)], (103) 

that is half a period out of phase with respect to Ug^x- When testing the hybrid 
method in two dimensions we use 4 particles per cell and a we interpolate the 
particle quantities to the grid using a Triangular-Shape-Cloud interpolation 
scheme [18]. 

The above initial conditions produce a sinusoidal wave with amplitude A prop- 
agating in the domain along the direction defined by the vector k. While we 
have experimented with various values for the parameters A, k and kq, below 
we present results for a few cases only, summarized in Table 3. In particular 
we consider a perturbation amplitude ^4 = 1.4 x 10-2 and adopt coupling co- 
efficients Ko = 1, 102, 10^ to explore the non-stiff and stiff regimes respectively 
and, for simplicity, we report only results from one dimensional tests for the 
two-fluid model and from two-dimensional tests for the hybrid model. 

In order to measure the rate at which the numerical solution converges, for 
each problem we carry out a set of 5 simulation runs employing N^eu = 
16, 32, 64, 128, 256 for a total range of 32. Note that the stiffness conditions do 
not change significantly as the grid is refined within the range of resolutions 
considered here. The convergence rate is measured using Richardson extrapo- 
lation. Given the numerical solution qr at a given resolution r we first estimate 
the error at a given point (i, j), as 

Sr;i,3 = <lr{hj) " 9r+l(^, j), (104) 



25 



where Qr+i is the solution at the next finer resolution, properly spatially aver- 
aged onto the coarser grid. We then take the n-norm of the error 



(105) 



where, Vij — Ax is the cell volume, and estimate the convergence rate as 

^" ^ ln(Ax./Aa;,) " 

For each studied case listed in Table 3, we produce a corresponding Table 4-7 
reporting the Li, L2 and Loo norms of the error and the corresponding conver- 
gence rates, Ri, R2 and Roo, as defined above. Note that the convergence rate 
of the dust component is evaluated in a way analogous to the fluid components, 
i.e. by considering the convergence of the error of the fluid representation of 
the particles. 



6.2.1 Two- fluid 

Cases A-E test the performance of the two-fluid scheme in which the dust is 
fully treated as a fluid using the scheme described in Sec. 3-4 and the parti- 
cle scheme is not used. The flve tests correspond to the following cases: (A) 
non-stiff, (B-D) stiff but with a range of dust-to-gas ratios and (E) mildly stiff 
case. In all tests the perturbation is along the x-axis, k = (1,0), although we 
have tested that the performance is the same in more than one dimension. The 
convergence rates for each case are reported in Tables (4)- (8), respectively, for 
the density and velocity of the gas and dust. Errors in the gas pressure are 
not reported but exhibit the same behavior. In the non-stiff case (A) all quan- 
tities converge with second order accuracy, except the dust density, which is 
somewhat expected since in the absence of drag a pressure-less fluid becomes 
singular. In the stiff cases (B-D), the scheme produces second order accu- 
rate results independent of the dust-to- gas ratio. In the intermediate regime, 
KoAt ~ 1, the convergence rate is between first and second order accurate 
as expected theoretically [19]. Thus, the two-fiuid algorithm for the dust-fiuid 
system is second order accurate in both the stiff and non-stiff regimes, and 
somewhat less accurate in between. 



6.2.2 Full Scheme 

Cases F-L test the performance of the full fluid-|-particle in which the dust 

is represented by a set of particles whose velocity and position are updated 
with the scheme described in Sec. 2. As in the previous section, the five tests 
correspond to the following cases: (F) non-stiff, (G-I) stiff but with a range 
of dust-to-gas ratios and (L) mildly stiff case. In all tests the perturbation is 
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Table 4 



Conver 


gence Rates: C 


ase : A = 


1.4 X 10-2, k = 


(1,0), 


Pd/Pg = 


1. 










Ncells 


Li 


Ri 




R2 


Loo 


Roo 


Li 


Ri 


L2 


R2 


Loo 


Roo 




density-gas 










x-vel-gas 










16 


6.4E-07 




1.4E-06 


- 


4.1E-06 




3.2E-06 




7.1E-06 




2.0E-05 




32 


1.6E-07 


2.0 


3.6E-07 


2.0 


l.OE-06 


2.0 


8.1E-07 


2.0 


1.8E-06 


2.0 


5.1E-06 


2.0 


64 


4.0E-08 


2.0 


8.8E-08 


2.0 


2.5E-07 


2.0 


2.0E-07 


2.0 


4.5E-07 


2.0 


1.3E-06 


2.0 


128 


9.9E-09 


2.0 


2.2E-08 


2.0 


6.3E-08 


2.0 


5.1E-08 


2.0 


l.lE-07 


2.0 


3.2E-07 


2.0 




density-dust 










x-vel-dust 










16 


1.3E-05 




2.9E-05 




8.2E-05 




4.6E-06 




l.OE-05 




2.9E-05 




32 


6.1E-06 


1.1 


1.4E-05 


1.1 


3.8E-05 


1.1 


1.2E-06 


2.0 


2.6E-06 


2.0 


7.2E-06 


2.0 


64 


2.9E-06 


1.0 


6.5E-06 


1.0 


1.9E-05 


1.0 


2.9E-07 


2.0 


6.4E-07 


2.0 


1.8E-06 


2.0 


128 


1.5E-06 


1.0 


3.2E-06 


1.0 


9.2E-06 


1.0 


7.2E-08 


2.0 


1.6E-07 


2.0 


4.5E-07 


2.0 



skewed with respect to the x-axis, k = (2, l)/^/5. The convergence rates for 
each case are reported in Tables (9)-(13), respectively, for all of the gas and 
dust variables. 

Inspection of the reported tables shows that, as expected, in the non-stiff 
regime (Table 9) the error drops with second order accuracy, this time even 
for the dust density. In the stiff regime and, unless the fluid mass dominates 
over the dust mass, the error on the fluid quantities drops approximately 
with flrst order accuracy, see Table 10-12. Since, as illustrated in the Sec. 6.1 
and 6.2.1, both the particle scheme and the two-fluid scheme retain their 
second order accuracy irrespective of the stiffness conditions, the worsening 
in accuracy is most likely due to the difficulty of coupling the gas and the 
dust fully self-consistently in the stiff regimes. However, the scheme is stable 
and convergent, though only flrst order, which would not have been trivially 
expected. Finally, in the mildly stiff regime, the convergence rate reduces to 
approximately flrst order, as shown in Tables 13 dust-to-gas ratio of 1. This 
is however, not unexpected because in this case accuracy of the two-fluid 
algorithm also drops. 
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Table 5 

Convergence Rates: Case : B = lA x 10~^, k = (1,0), Pd/Pg = 1- 



Nceiis Li Ri L2 R2 -^00 Roo Li Ri L2 R2 R(x, 

density-gas x-vel-gas 

16 3.2E-06 - 7.1E-06 - 2.1E-05 - 9.3E-07 - 2.3E-06 - 9.6E-06 - 

32 8.0E-07 2.0 1.8E-06 2.0 5.2E-06 2.0 2.1E-07 2.1 4.9E-07 2.2 2.3E-06 2.1 

64 2.0E-07 2.0 4.4E-07 2.0 1.3E-06 2.1 5.7E-08 1.9 1.3E-07 1.9 4.5E-07 2.3 

128 4.9E-08 2.0 l.lE-07 2.0 3.1E-07 2.0 1.6E-08 1.8 3.6E-08 1.8 1.2E-07 1.9 



density-dust x-vel-dust 

16 3.7E-05 - 8.9E-05 - 4.0E-04 - 1.3E-06 - 2.9E-06 - 9.1E-06 - 

32 l.OE-05 1.9 2.5E-05 1.8 1.3E-04 1.6 3.4E-07 1.9 7.6E-07 1.9 3.4E-06 1.4 

64 2.6E-06 2.0 6.0E-06 2.1 2.0E-05 2.7 7.3E-08 2.2 1.7E-07 2.2 7.3E-07 2.2 

128 7.7E-07 1.8 1.7E-06 1.8 4.8E-06 2.1 1.8E-08 2.0 4.1E-08 2.0 1.9E-07 1.9 

Table 6 

Convergence Rates: Case : C = 1.4 x 10"^, k = (1,0), Pa/Pg = 10~^. 

Ncells Li Ri L2 R2 -^00 -Roo Li Ri L2 R2 Ltx, Roo 



density-gas x-vel-gas 

16 4.2E-07 - 9.5E-07 - 2.7E-06 - 2.5E-06 ~ 5.7E-06 ~ 1.6E-05 

32 l.lE-07 2.0 2.4E-07 2.0 6.8E-07 2.0 6.4E-07 2.0 1.4E-06 2.0 4.0E-06 2.0 

64 2.7E-08 2.0 6.0E-08 2.0 1.7E-07 2.0 1.6E-07 2.0 3.6E-07 2.0 l.OE-06 2.0 

128 6.8E-09 2.0 1.5E-08 2.0 4.3E-08 2.0 4.0E-08 2.0 9.0E-08 2.0 2.5E-07 2.0 



density-dust x-vel-dust 

16 2.5E-08 - 5.6E-08 - 1.6E-07 - 1.7E-06 - 3.7E-06 - l.lE-05 - 

32 4.8E-09 2.4 l.lE-08 2.4 3.0E-08 2.4 4.2E-07 2.0 9.3E-07 2.0 2.6E-06 2.0 

64 3.6E-10 3.7 8.0E-10 3.7 2.4E-09 3.7 l.lE-07 2.0 2.3E-07 2.0 6.7E-07 2.0 

128 3.5E-10 0.1 7.7E-10 0.1 2.2E-09 0.1 2.7E-08 2.0 5.9E-08 2.0 1.7E-07 2.0 
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Table 7 

Convergence Rates: Case : D = 1.4 x 10~^, k = (1,0), pa/ pg = 10^. 



Nceiis Li Ri L2 R2 -^00 Roo Li Ri L2 R2 R(x, 

density-gas x-vel-gas 

16 3.1E-07 - 6.9E-07 - 1.9E-06 - 5.5E-06 - 1.2E-05 - 3.4E-05 - 

32 7.5E-08 2.0 1.7E-07 2.0 4.7E-07 2.0 1.4E-06 2.0 3.1E-06 2.0 8.7E-06 2.0 

64 1.8E-08 2.1 4.0E-08 2.1 l.lE-07 2.1 3.5E-07 2.0 7.8E-07 2.0 2.2E-06 2.0 

128 4.1E-09 2.1 9.2E-09 2.1 2.6E-08 2.1 8.8E-08 2.0 1.9E-07 2.0 5.5E-07 2.0 



density-dust x-vel-dust 

16 1.4E-06 - 3.0E-06 - 8.5E-06 - 5.5E-06 - 1.2E-05 - 3.4E-05 - 

32 3.3E-07 2.0 7.4E-07 2.0 2.1E-06 2.0 1.4E-06 2.0 3.1E-06 2.0 8.7E-06 2.0 

64 7.8E-08 2.1 1.7E-07 2.1 4.9E-07 2.1 3.5E-07 2.0 7.8E-07 2.0 2.2E-06 2.0 

128 1.7E-08 2.2 3.8E-08 2.2 l.lE-07 2.2 8.8E-08 2.0 1.9E-07 2.0 5.5E-07 2.0 

Table 8 

Convergence Rates: Case : E = lA x 10~^, k = (1,0), Pd/Pg = 1- 

Ncells Li Ri L2 R2 -^00 -Roo Li Ri L2 R2 Ltx, Roo 



density-gas x-vel-gas 

16 3.0E-06 - 6.8E-06 - 1.9E-05 - 2.6E-06 - 5.9E-06 - 1.7E-05 

32 1.5E-06 1.0 3.3E-06 1.0 9.4E-06 1.0 l.lE-06 1.3 2.4E-06 1.3 6.7E-06 1.3 

64 7.4E-07 1.0 1.6E-06 1.0 4.6E-06 1.0 4.0E-07 1.4 8.8E-07 1.4 2.5E-06 1.4 

128 2.3E-07 1.7 5.1E-07 1.7 1.5E-06 1.7 1.2E-07 1.8 2.6E-07 1.8 7.4E-07 1.8 



density-dust x-vel-dust 

16 3.3E-05 - 7.3E-05 - 2.0E-04 - 2.7E-06 - 6.0E-06 - 1.7E-05 - 

32 9.0E-06 1.9 2.0E-05 1.9 5.6E-05 1.9 l.OE-06 1.4 2.3E-06 1.4 6.5E-06 1.4 

64 3.0E-06 1.6 6.7E-06 1.6 1.9E-05 1.6 3.6E-07 1.5 8.1E-07 1.5 2.3E-06 1.5 

128 8.7E-07 1.8 1.9E-06 1.8 5.5E-06 1.8 l.OE-07 1.8 2.3E-07 1.8 6.6E-07 1.8 
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Table 9 

Convergence Rates: Case : F = 1.4 x 10~^, k = (2, l)/-\/5, Pd/ Pg = 1- 



Ncells Li Ri L2 R2 L(x, Roo Li Ri L2 R2 Ltx, Roo 

density-gas x-vel-gas 

16 2.1E-05 - 2.4E-05 - 3.3E-05 - 5.4E-05 - 6.0E-05 - 8.5E-05 - 

32 6.0E-06 1.8 6.6E-06 1.8 9.4E-06 1.8 1.5E-05 1.9 1.6E-05 1.9 2.3E-05 1.9 

64 1.6E-06 1.9 1.7E-06 1.9 2.5E-06 1.9 3.8E-06 2.0 4.2E-06 2.0 6.0E-06 2.0 

128 4.0E-07 2.0 4.4E-07 2.0 6.3E-07 2.0 9.6E-07 2.0 l.lE-06 2.0 1.5E-06 2.0 



y-vel-gas pressure 

16 1.7E-05 - 1.8E-05 - 2.6E-05 - 3.0E-05 - 3.3E-05 - 4.7E-05 - 

32 4.8E-06 1.8 5.4E-06 1.8 7.6E-06 1.8 8.4E-06 1.8 9.3E-06 1.8 1.3E-05 1.8 

64 1.3E-06 1.9 1.4E-06 1.9 2.0E-06 1.9 2.2E-06 1.9 2.4E-06 1.9 3.4E-06 1.9 

128 3.3E-07 2.0 3.6E-07 2.0 5.1E-07 2.0 5.6E-07 2.0 6.2E-07 2.0 8.8E-07 2.0 



density-dust x-vel-dust 

16 2.3E-04 - 2.5E-04 - 3.6E-04 - 7.8E-05 - 8.6E-05 - 1.2E-04 

32 4.9E-05 2.2 5.5E-05 2.2 8.0E-05 2.2 2.0E-05 1.9 2.2E-05 1.9 3.2E-05 1.9 

64 1.6E-05 1.7 1.7E-05 1.7 2.5E-05 1.7 5.1E-06 2.0 5.7E-06 2.0 8.0E-06 2.0 

128 3.8E-06 2.1 4.3E-06 2.0 9.4E-06 1.4 1.3E-06 2.0 1.4E-06 2.0 2.0E-06 2.0 



y-vel-dust 

16 1.3E-04 - 1.5E-04 - 2.1E-04 

32 3.5E-05 1.9 3.9E-05 1.9 5.5E-05 1.9 

64 8.8E-06 2.0 9.8E-06 2.0 1.4E-05 2.0 

128 2.2E-06 2.0 2.5E-06 2.0 3.5E-06 2.0 



6.3 Streaming Instability 

In this last section we carry out the test proposed in [28] (hereafter, YJ07), 
which consists of following the growth of an initial perturbation unstable to the 
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Table 10 

Convergence Rates: Case : G = 1.4 x 10~^, k = (2, l)/-\/5, Pd/ Pg = 1- 



Nceiis Li Ri L2 R2 L(X) Roo Li R\ L2 R2 R, 



'OO 



density-gas x-vel-gas 

16 5.8E-05 - 6.4E-05 - 9.7E-05 - 5.0E-05 - 5.6E-05 - 8.7E-05 - 

32 1.3E-05 2.2 1.5E-05 2.1 2.4E-05 2.0 1.6E-05 1.6 1.8E-05 1.6 2.8E-05 1.6 

64 3.7E-06 1.8 4.1E-06 1.9 6.3E-06 1.9 5.3E-06 1.6 6.0E-06 1.6 8.7E-06 1.7 

128 1.5E-06 1.3 1.6E-06 1.3 2.6E-06 1.3 2.0E-06 1.4 2.3E-06 1.4 3.5E-06 1.3 

y-vel-gas pressure 

16 2.5E-05 - 3.0E-05 - 7.0E-05 - 8.2E-05 - 9.0E-05 - 1.4E-04 - 

32 3.2E-05 -0.4 3.6E-05 -0.2 5.6E-05 0.3 1.8E-05 2.2 2.1E-05 2.1 3.3E-05 2.0 

64 2.5E-05 0.4 2.8E-05 0.4 3.9E-05 0.5 5.2E-06 1.8 5.8E-06 1.9 8.8E-06 1.9 

128 1.5E-05 0.7 1.6E-05 0.7 2.3E-05 0.8 2.1E-06 1.3 2.3E-06 1.3 3.7E-06 1.3 

density-dust x-vel-dust 

16 9.6E-05 - l.lE-04 - 2.0E-04 - 4.9E-05 - 5.4E-05 - 8.1E-05 

32 3.8E-05 1.3 4.4E-05 1.3 7.8E-05 1.3 1.8E-05 1.4 2.0E-05 1.4 2.9E-05 1.5 

64 1.7E-05 1.1 1.9E-05 1.2 3.1E-05 1.3 6.1E-06 1.6 6.8E-06 1.6 9.8E-06 1.6 

128 8.0E-06 1.1 9.1E-06 1.1 1.5E-05 1.1 2.2E-06 1.4 2.5E-06 1.4 3.8E-06 1.4 

y-vel-dust 

16 1.8E-04 - 2.0E-04 - 2.9E-04 
32 8.1E-05 1.2 9.0E-05 1.2 1.3E-04 1.2 
64 3.8E-05 1.1 4.2E-05 1.1 5.9E-05 1.1 
128 1.8E-05 1.1 2.0E-05 1.1 2.8E-05 1.1 

streaming instability [29]. The test setup, while considerably simplified with 

respect to a Keplerian disk, preserves the key dynamical features leading to the 
streaming instability in a protoplanetary disk [14,15] and, with the available 
analytic solution (for the linear regime), provides a clean test for a simula- 
tion code. The test consists of simulating an axisymmetric shearing sheet in 
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Table 11 

Convergence Rates: Case : H = lA x 10~^, k = (2, l)/-/5, Pd/Pg = 10~^. 



Nceiis Li Ri L2 R2 L(x, Roo Li Ri L2 R2 Ltx, R, 



'OO 



density-gas x-vel-gas 

16 1.5E-05 - 1.7E-05 - 2.4E-05 - 4.3E-05 - 4.8E-05 - 6.7E-05 - 

32 4.5E-06 1.8 5.0E-06 1.8 7.1E-06 1.8 1.2E-05 1.9 1.3E-05 1.9 1.9E-05 1.8 

64 1.2E-06 1.9 1.3E-06 1.9 1.9E-06 1.9 3.0E-06 2.0 3.4E-06 2.0 4.8E-06 2.0 

128 3.0E-07 2.0 3.3E-07 2.0 4.8E-07 2.0 7.5E-07 2.0 8.4E-07 2.0 1.2E-06 2.0 

y-vel-gas pressure 

16 l.lE-06 - 1.2E-06 - 1.8E-06 - 2.1E-05 - 2.4E-05 - 3.4E-05 - 

32 5.4E-07 1.0 6.0E-07 1.1 8.6E-07 1.1 6.3E-06 1.8 7.0E-06 1.8 9.9E-06 1.8 

64 2.0E-07 1.4 2.3E-07 1.4 3.3E-07 1.4 1.7E-06 1.9 1.8E-06 1.9 2.6E-06 1.9 

128 6.4E-08 1.7 7.1E-08 1.7 l.OE-07 1.7 4.2E-07 2.0 4.7E-07 2.0 6.7E-07 2.0 

density-dust x-vel-dust 

16 4.1E-07 - 4.6E-07 - 6.5E-07 - 1.8E-04 - 2.0E-04 - 2.8E-04 

32 1.4E-07 1.6 1.5E-07 1.6 2.1E-07 1.6 4.9E-05 1.9 5.4E-05 1.9 7.6E-05 1.9 

64 5.2E-08 1.4 5.8E-08 1.4 8.2E-08 1.4 1.2E-05 2.0 1.4E-05 2.0 1.9E-05 2.0 

128 2.0E-08 1.4 2.2E-08 1.4 3.4E-08 1.3 3.1E-06 2.0 3.5E-06 2.0 4.9E-06 2.0 

y-vel-dust 

16 3.4E-05 - 3.7E-05 - 5.2E-05 
32 9.9E-06 1.8 l.lE-05 1.8 1.6E-05 1.8 
64 2.6E-06 1.9 2.9E-06 1.9 4.1E-06 1.9 
128 6.7E-07 2.0 7.4E-07 2.0 l.OE-06 2.0 

which dust particles drift with respect to the gas due to the radial pressure 
gradient affecting the gas but not the dust. Both vertical structure (along the 
y-direction) and self-gravity are ignored. After neglecting the vertical deriva- 
tives (along the y-axis), the governing equations for the gas component in a 
shearing sheet can be written as in Eq. (1) provided that we replace fa in the 
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Table 12 

Convergence Rates: Case : / = 1.4 x 10~^, k = (2, l)/-\/5, Pdl Pg = 10^. 



Ncells 



Ri L2 



R2 



R^ 



L2 R2 



Lr 



Rr 



density-gas x-vel-gas 

16 6.5E-06 - 7.2E-06 - l.OE-05 - 1.7E-04 - 1.9E-04 - 2.6E-04 - 

32 2.3E-06 1.5 2.5E-06 1.5 3.5E-06 1.5 1.4E-04 0.3 1.6E-04 0.3 2.2E-04 0.2 

64 l.OE-06 1.2 l.lE-06 1.2 1.6E-06 1.2 8.9E-05 0.7 9.9E-05 0.7 1.4E-04 0.7 

128 5.3E-07 0.9 5.9E-07 0.9 8.3E-07 0.9 5.4E-05 0.7 6.0E-05 0.7 8.4E-05 0.7 

y-vel-gas pressure 

16 2.3E-04 - 2.6E-04 - 3.6E-04 - 9.1E-06 - l.OE-05 - 1.4E-05 - 

32 2.0E-04 0.2 2.2E-04 0.2 3.1E-04 0.2 3.2E-06 1.5 3.5E-06 1.5 5.0E-06 1.5 

64 1.2E-04 0.7 1.4E-04 0.7 2.0E-04 0.6 1.4E-06 1.2 1.5E-06 1.2 2.2E-06 1.2 

128 7.5E-05 0.7 8.3E-05 0.7 1.2E-04 0.7 7.4E-07 0.9 8.2E-07 0.9 1.2E-06 0.9 

density-dust x-vel-dust 

16 2.2E-05 - 2.4E-05 - 3.3E-05 - 3.8E-04 - 4.1E-04 - 5.8E-04 

32 9.1E-06 1.2 l.OE-05 1.2 1.4E-05 1.2 2.0E-04 0.9 2.3E-04 0.9 3.2E-04 0.9 

64 4.0E-06 1.2 4.4E-06 1.2 6.2E-06 1.2 l.lE-04 0.9 1.2E-04 0.9 1.7E-04 0.9 

128 1.9E-06 1.1 2.1E-06 1.1 3.0E-06 1.1 5.8E-05 0.9 6.5E-05 0.9 9.2E-05 0.9 



y-vel-dust 

16 5.2E-04 ~ 5.8E-04 - 8.1E-04 

32 2.9E-04 0.9 3.2E-04 0.9 4.5E-04 0.8 

64 1.5E-04 0.9 1.7E-04 0.9 2.4E-04 0.9 

128 8.2E-05 0.9 9.1E-05 0.9 1.3E-04 0.9 



source term in (2) with 



2 "'9,x 



(107) 
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Table 13 

Convergence Rates: Case : L = 1.4 x 10"^ k = (2, 1)/V5, pd/ pg = 1- 



Nceiis Li Ri L2 R2 L(X) Roo Li Ri L2 R2 Ltx, Roo 

density-gas x-vel-gas 

16 3.2E-05 - 3.5E-05 - 4.9E-05 - 4.3E-05 - 4.8E-05 - 6.7E-05 

32 8.1E-06 2.0 9.0E-06 2.0 1.3E-05 1.9 4.4E-05 -0.0 4.9E-05 -0.0 6.9E-05 -0.0 

64 1.6E-05 -1.0 1.8E-05 -1.0 2.5E-05 -1.0 4.3E-05 0.0 4.8E-05 0.0 6.8E-05 0.0 

128 8.8E-06 0.9 9.7E-06 0.9 1.4E-05 0.9 1.9E-05 1.2 2.1E-05 1.2 3.0E-05 1.2 

y-vel-gas pressure 

16 2.6E-05 - 2.9E-05 - 4.1E-05 - 4.4E-05 - 4.9E-05 - 6.9E-05 

32 8.1E-05 -1.6 9.0E-05 -1.6 1.3E-04 -1.6 l.lE-05 2.0 1.3E-05 2.0 1.8E-05 2.0 

64 5.9E-05 0.5 6.6E-05 0.5 9.3E-05 0.5 2.2E-05 -1.0 2.5E-05 -1.0 3.5E-05 -1.0 

128 2.4E-05 1.3 2.7E-05 1.3 3.8E-05 1.3 1.2E-05 0.9 1.4E-05 0.9 1.9E-05 0.9 

density-dust x-vel-dust 

16 1.3E-04 - 1.4E-04 - 2.0E-04 - 5.2E-05 - 5.7E-05 - 8.0E-05 

32 1.6E-04 -0.3 1.7E-04 -0.3 2.5E-04 -0.3 4.1E-05 0.3 4.6E-05 0.3 6.5E-05 0.3 

64 l.lE-04 0.6 1.2E-04 0.5 1.7E-04 0.5 4.4E-05 -0.1 4.8E-05 -0.1 6.9E-05 -0.1 

128 4.4E-05 1.3 4.9E-05 1.3 7.0E-05 1.3 2.0E-05 1.2 2.2E-05 1.2 3.1E-05 1.2 

y-vel-dust 

16 1.8E-04 - 2.0E-04 - 2.8E-04 

32 1.2E-04 0.5 1.4E-04 0.5 1.9E-04 0.5 

64 6.9E-05 0.8 7.7E-05 0.8 l.lE-04 0.8 

128 2.7E-05 1.4 3.0E-05 1.4 4.2E-05 1.4 



where Q is the angular velocity of the shearing sheet, V}- — is the azimuthal 
velocity at radial distance r, and 7] is a, dimensionless parameter expressing the 
strength of the radial pressure gradient with respect to the centrifugal force. 
The source term for the dust fluid is modified in an analogous way, except 
that in this case we set 77 = 0. Finally, the particle equations of motion are 
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also modified, namely Eq. (4) now reads 



dt 



CO 



^ 2 0^ 

-I 

y Oy 



(108) 



In principle the y-component of the particles position should also be modified, 
but this is irrelevant due to the azimuthal symmetry of the system. 
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Table 14 
Run Set 



case 




k 


Pd/pg 


s/Vt 


linA 


0.357 


(30,0,30) 


3 


0.4190204 


linB 


0.357 


(6,0,6) 


0.2 


0.0154764 



In order to perform the test, the gas and dust velocities are initialized to the 
equilibrium solution given by Nakagawa et al. [22] (or Eq. 7 in YJ05). Gas 
density and pressure are set according to 

p, = 1.4, P=-^, (109) 

where Xp = rjVk/c = 0.05 as in YJ05, and pa is set according to a specified 
dust-to-gas ratio (sec below). Note that the parameters Q, r and rj need not 
be specified as long as one expresses time, distances and velocities in units 
of rjr, and rjflr = r]Vk, respectively. Finally, a perturbation unstable to 
streaming instability, is added to the initial equilibrium values. YJ07 provide 
two sets of the eigenmodes for the perturbation amplitudes of each variable, 
corresponding to what they refer to as linA and linB cases. Some basic fea- 
tures of these cases, including the wave vector, the dust-to-gas ratio and the 
normalized growth rate {s/Q), are summarized in Table 14, full details are 
provided in YJ07. The perturbations are in the radial and vertical directions 
(x,z) but not in the azimuthal direction (y). The system is evolved using an 
isothermal equation of state. Note that these tests are basically in the non-stiff 
regime, as kg^d^t < 0.1. 

In Fig. 3 we report the simulated growth rate of the instabihty, for both gas 
and dust variables, as a function of the resolution expressed in number of cells 
per perturbation wavelength. The top and bottom panels correspond to case 
linA and linB, respectively, and open and filled dots refer to the two-fiuid 
and hybrid algorithm, respectively. With respect to the hybrid algorithm, we 
always employ 16 particles per cell and a TSC interpolation scheme. As for the 
linA case, the test indicates that for both the two-fluid and hybrid algorithms, 
at least 32 cells per wavelength are required in order to capture the instability 
growth for all variables and convergence to the analytic solution is basically 
achieved with 64 cells per wavelength. The linB case is much more challenging 
and now at least 128 cells per wavelength are necessary in order to capture the 
instabihty. Note so that because the gas-dust coupling is non-stiff, we do not 
expect any particular advantage of our scheme with respect to other second 
order schemes. In fact, our results are comparable to those of Balsara et al. [1], 
who use an approach based on a second order Godunov's method, and less 
accurate then those in YJ05, who instead use a sixth order spectral method. 
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7 Discussion &c Summary 



We have presented a stable and convergent method for studying a system of 
gas and dust coupled through viscous drag in both non-stiff and stiff regimes. 
Our approach consists of updating the fluid quantities using a two fluid model 
and then using the updated fluid solution to advance the individual parti- 
cle solutions with a self- consistent time evolution of the gas velocity in the 
estimate of the drag force. 

In our derivation of the two fluid method, we first obtain a fluid description 
of the dust component using a Particle-Mesh method, and then we study 
the modified gas-dust hyperbolic system following the approach in Miniati & 
Colella [19]. Based on this analysis we formulate a predictor step providing 
first order accurate reconstruction of the time-averaged state variables at cell 
interfaces, whence a second order accurate estimates of the conservative fluxes 
can be obtained. Finally, for the time-discretization for the source terms we use 
a single-step, second-order accurate scheme derived from the a— QSS method 
proposed by Mott, Oran and van Leer ([21]). This completes the description 
of our two fluid method. 

The fluid description of the dust component (Eq. 23, 24) assumes the simplest 
type of closure which neglects dispersion velocity terms in the momentum 
equation. However, since the particle distribution is known, different and more 
suitable closures can be constructed according to the specifics of the aimed 
application. This is relevant when the dust-gas coupling is stiff and the dust 
backreaction is important, e.g. high dust density and large dust grain size [13], 
because in this case the gas tends to follow the dust and the anisotropic particle 
motions may not be quickly damped. 

In order to advance the individual particle solutions we use the fluid solution 
to determine the time dependence of the gas velocity entering the drag terms. 
This allows us to derive a particle integration scheme, also based on the a-QSS 
method, which is second order accurate regardless of the strength of dust-gas 
coupling. Remarkably, the particle method that we employ contains no explicit 
term arising from the stiff coupling of the particle component, in the sense that 
each particle motion is integrated individually, though it effectively depends 
on the other particles solutions, and the scheme is essentially explicit in time, 
as it only involves the particle solution at time t = nAt. Note that when the 
dust backreaction on the gas is not important (low dust density) , the particle 
scheme presented in Sec. 2 can be used together with an unmodified Godunov's 
method for the gas component. Similarly, a simpler scheme can be obtained in 
the limit Kg — )■ 0, as the gas description of the semi-implicit two-fluid scheme 
in Sec. 3 reduces to an ordinary explicit Godunov's method, while the stiff 
solver remains unchanged for the dust component. 
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A set of benchmark problems show that our method is stable and convergent. 
In particular, the two-fluid approach is second order accurate both in the non- 
stifi^ and stiff regimes, and drops to first order in the intermediate regime as 
expected theoretically [19]. The hybrid scheme, on the other hand, is second 
order only in the non-stiff regime. Since, as illustrated in the Sec. 6.1 and 6.2.1, 
both the particle scheme and the two-fluid scheme retain their second order 
accuracy irrespective of the stiffness conditions, the drop in accuracy in the 
hybrid approach is most likely due to the difficulty of coupling the gas and 
the dust fully self-consistently in the stiff regimes. At any rate, the scheme 
remains stable and first order convergent even in the stiff regime. 

We have also tested our code against the streaming instability problem pre- 
sented in YJ05. This is a clean test and a very relevant one for protoplanetary 
disk applications. However, due to the specifics of the test set up, the gas-dust 
coupling is non-stiff, so that we do not expect any particular advantage of our 
scheme with respect to other second order schemes. In fact, our results are 
comparable to those of Balsara et al. [1], who use an approach based on a 
second order Godunov's method, and less accurate then those in YJ05, who 
instead use a sixth order spectral method. 

Finally, although the present analysis focuses on the Epstein regime for the 
functional form of drag force, in principle extension to the case of Stokes regime 
should also be possible. This would require a modification of both the Godunov 
predictor step and the semi-implicit scheme in Sec. 3.2. In addition, the present 
scheme can also be extended to the case of dust particles of multiple sizes. A 
complication arises, however, due to the fact that one cannot straightforwardly 
define a single effective dust component to which the gas is coupled via the drag 
force. Such complication seems inevitable and of general character, i.e. it is not 
restricted to the approach presented in this paper, when the coupling is stiff 
and the dust backreaction is dynamically important. So, for a limited number 
of dust components of different grain size the present method can be modified 
by defining an extended system in which each dust component is represented 
separately. This approach, however, becomes cumbersome and expensive when 
the number of dust species is large, and a method for constructing a single 
effective dust component should be investigated in such cases. 
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A Derivation of the a-QSS Based Pcirticle Scheme 



We use the alpha-QSS method to integrate the following equations, describing 
the trajectories in phase-space of dust particles, 



dxd 

-^-Vd, (A.l) 



'd 



"dt 



-^^d {vd - u,) - V(/., (A.2) 



from time t — i", when x^ii) — x^, Vd{t) — v^, tot — = t"-+At. Following 
Eq. (11) we first integrate Eq. (A.2) assuming po = K'd ^"<i Qo — i^^u^ — V0, 
which yields 

Vd{t) = v2 + « - ^^^)(1 - e-'^S*) _ t V0. (A.3) 

Note that we have modified the gravitational acceleration term because, as 
discussed in Sec. 3.2, gravitational acceleration is not attenuated by drag as 
both gas and dust are equally accelerated by it. We then obtain a first order 
accurate expression of the particle position by time integration of the above 
solution. 



Xd{At)=x'l + dTVdir)^ 

= + Atvl + At{u- - v2) [1 - ^^^^ J - — V0, (A.4) 

which is equivalent to Eq. (17). For the final corrector step we need to estimate 
the equivalent of and p appearing in (12). The expressions for 5* and and p 

in Eq. (13) and (14), respectively, can be obtained by assuming a linear time 
dependence for these parameters. Instead for the gas velocity (which play the 
equivalent role of q in our case) we use the ansatz in Eq. (18) with 

A^^ = v^+' - vl 

On the other hand, for the drag coefficients (which play the equivalent role of 
p in our case) we use the following averages consistent with the prescription 
in Eq. (14), 

'^s = ^[i^s{Xd) + i^sixd)], s = d,g. 

Using Eq. (18) and the above expressions for Avy and k^, Kg, we can then 
integrate Eq. (A.2) to obtain a time dependent solution of the dust particle 
velocity. Its evaluation at time t"'~^^ gives the solution in Eq. (20) and its 
integration over a timestep At gives the solution (19). 
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